library(nnet)
library(tidyverse)
library(marginaleffects)
library(ggpubr)
library(stargazer)
library(haven)
library(psych)
rm(list=ls())

df<-read_sav("fact_opinion_diff.sav")


# Matrix of response types
df$Y<-matrix(c(df$accurate_response,df$partisan_error,df$unbiased_error),ncol=3,dimnames = list(c(),c("accurate_response","partisan_error","unbiased_error")))

# Model for first figure
step1<-multinom(Y ~ polknow+analytical+educ+currevent+age+repdum+demdum,Hess=TRUE, data = df)


# Response predictions for new data

# Civics knowledge
pk<-predictions(
  step1,
  newdata = datagrid(polknow=c(0:6),repdum=c(0),demdum=c(0)),
  vcov = TRUE,
  conf_level = 0.95)

# cognitive  ability
an<-predictions(
  step1,
  newdata = datagrid(analytical=c(0:8),repdum=c(0),demdum=c(0)),
  vcov = TRUE,
  conf_level = 0.95)

# Education
ed<-predictions(
  step1,
  newdata = datagrid(educ=c(0:5),repdum=c(0),demdum=c(0)),
  vcov = TRUE,
  conf_level = 0.95)

#Current events
ce<-predictions(
  step1,
  newdata = datagrid(currevent=c(1:10),repdum=c(0),demdum=c(0)),
  vcov = TRUE,
  conf_level = 0.95)

pk$Outcomes<-pk$group
an$Outcomes<-an$group
ed$Outcomes<-ed$group
ce$Outcomes<-ce$group

col<-c("green4","red","blue")
pk1<-pk%>%
ggplot(aes(x = polknow, y = estimate)) +   
  geom_line(aes(colour = Outcomes), linewidth = .5) +
  xlab("Civics knowledge") +
  ylab("")+
  ylim(0,1)+
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = Outcomes), 
              alpha = 0.2) + 
  scale_fill_manual(values=col)+
  scale_color_manual(values=col)+
  theme(text = element_text(size = 9), 
    panel.grid.major.y = element_line(color="grey",
                                          linewidth = .5,
                                          linetype=2),
    axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
    legend.background = element_rect(fill="white"),
        plot.background = element_rect(fill="white"),
        panel.background = element_rect(fill="white"),
        panel.border = element_rect(color="grey",linewidth=1,fill="transparent"))

an1<-an%>%
  ggplot(aes(x = analytical, y = estimate)) +   
  geom_line(aes(colour = Outcomes), linewidth = .5) +
  xlab("Cognitive ability") +
  ylab("")+
  ylim(0,1)+
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = Outcomes), 
              alpha = 0.2) + 
  scale_fill_manual(values=col)+
  scale_color_manual(values=col)+
  theme(text = element_text(size = 9),   
    panel.grid.major.y = element_line(color="grey",
                                          linewidth = .5,
                                          linetype=2),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        legend.background = element_rect(fill="white"),
        plot.background = element_rect(fill="white"),
        panel.background = element_rect(fill="white"),
        panel.border = element_rect(color="grey",linewidth=1,fill="transparent"))



ed1<-ed%>%
  ggplot(aes(x = educ, y = estimate)) +   
  geom_line(aes(colour = Outcomes), linewidth = .5) +
  xlab("Education") +
  ylab("")+
  ylim(0,1)+
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = Outcomes), 
              alpha = 0.2) + 
  scale_fill_manual(values=col)+
  scale_color_manual(values=col)+
  theme(text = element_text(size = 9), 
    panel.grid.major.y = element_line(color="grey",
                                          linewidth = .5,
                                          linetype=2),
        axis.text.y=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        legend.background = element_rect(fill="white"),
        plot.background = element_rect(fill="white"),
        panel.background = element_rect(fill="white"),
        panel.border = element_rect(color="grey",linewidth=1,fill="transparent"))


ce1<-ce%>%
  ggplot(aes(x = currevent, y = estimate)) +   
  geom_line(aes(colour = Outcomes), linewidth = .5) +
  xlab("Current events knowledge") +
  ylab("")+
  ylim(0,1)+
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = Outcomes), 
              alpha = 0.2) + 
  scale_fill_manual(values=col)+
  scale_color_manual(values=col)+
  theme(text = element_text(size = 9), 
    panel.grid.major.y = element_line(color="grey",
                                          linewidth = .5,
                                          linetype=2),
        axis.text.y=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        legend.background = element_rect(fill="white"),
        plot.background = element_rect(fill="white"),
        panel.background = element_rect(fill="white"),
        panel.border = element_rect(color="grey",linewidth=1,fill="transparent"))



################### 
# Creates first figure in paper
figure1<-ggarrange(pk1,ce1,an1,ed1,ncol=2, nrow = 2,legend="bottom",common.legend = TRUE)
annotate_figure(figure1,top = text_grob("Figure 1: Partisan Bias", face = "italic", size = 12))
###################

step2<-multinom(Y ~ dem_bias+rep_bias+polknow+analytical+educ+currevent+age+repdum+demdum, data = df,Hess = TRUE)



# Rep-rep bias
n.step1r <- expand.grid(0,0:100,mean(df$polknow,na.rm =TRUE),mean(df$analytical,na.rm =TRUE),mean(df$educ,na.rm =TRUE),mean(df$currevent,na.rm =TRUE),
                        mean(df$age,na.rm =TRUE),1,0)
colnames(n.step1r)<-c("dem_bias","rep_bias","polknow","analytical","educ","currevent","age","repdum","demdum")

# dem-dem bias
n.step1d <- expand.grid(0:100,0,mean(df$polknow,na.rm =TRUE),mean(df$analytical,na.rm =TRUE),mean(df$educ,na.rm =TRUE),mean(df$currevent,na.rm =TRUE),
                        mean(df$age,na.rm =TRUE),0,1)
colnames(n.step1d)<-c("dem_bias","rep_bias","polknow","analytical","educ","currevent","age","repdum","demdum")

dbias<-predictions(
  step2,
  newdata=n.step1d,
  vcov = TRUE,
  conf_level = 0.95)
rbias<-predictions(
  step2,
  newdata = n.step1r,
  vcov = TRUE,
  conf_level = 0.95)
rbias$Outcomes<-rbias$group
dbias$Outcomes<-dbias$group
col<-c("green4","red","blue")
db<-dbias%>%
  ggplot(aes(x = dem_bias, y = estimate)) +   
  geom_line(aes(colour = Outcomes), linewidth = .5) +
  xlab("Pro-Democrat affective polarization") +
  ylab("")+
  ylim(0,1)+
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = Outcomes), 
              alpha = 0.2) + 
  scale_fill_manual(values=col)+
  scale_color_manual(values=col)+
  theme(text = element_text(size = 9), 
        panel.grid.major.y = element_line(color="grey",
                                          linewidth = .5,
                                          linetype=2),
        #axis.text.x=element_blank(),
        #axis.ticks.x=element_blank(),
        plot.background = element_rect(fill="white"),
        panel.background = element_rect(fill="white"),
        panel.border = element_rect(color="grey",linewidth=1,fill="transparent"))

rb<-rbias%>%
  ggplot(aes(x = rep_bias, y = estimate)) +   
  geom_line(aes(colour = Outcomes), linewidth = .5) +
  xlab("Pro-Republican affective polarization") +
  ylab("")+
  ylim(0,1)+
  geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = Outcomes), 
              alpha = 0.2) + 
  scale_fill_manual(values=col)+
  scale_color_manual(values=col)+
  theme(text = element_text(size = 9), 
        panel.grid.major.y = element_line(color="grey",
                                          linewidth = .5,
                                          linetype=2),
        axis.text.y=element_blank(),
        #axis.text.x=element_blank(),
        #axis.ticks.x=element_blank(),
        plot.background = element_rect(fill="white"),
        panel.background = element_rect(fill="white"),
        panel.border = element_rect(color="grey",linewidth=1,fill="transparent"))


###############
# figure 2 in paper 
figure2<-ggarrange(db,rb,ncol=2, nrow = 1,legend="bottom",common.legend = TRUE)
###############

################################################################################
# APPENDIX 
################################################################################

###############
# figure 1
df$Fact1<-ifelse(df$FACT1==1,1,0) #FACT1  # Fact 
df$Fact2<-ifelse(df$FACT2==2,1,0)#FACT2  # opinion 
df$Fact3<-ifelse(df$FACT3==1,1,0)#FACT3  # Fact 
df$Fact4<-ifelse(df$FACT4==2,1,0)#FACT4  # opinion 
df$Fact5<-ifelse(df$FACT5==2,1,0)#FACT5  # opinion 
df$Fact6<-ifelse(df$FACT6==1,1,0)#FACT6  # Fact 
df$Fact7<-ifelse(df$FACT7==1,1,0)#FACT7  # Fact 
df$Fact8<-ifelse(df$FACT8==2,1,0)#FACT8  # opinion 
df$Fact9<-ifelse(df$FACT9==1,1,0)#FACT9  # Fact 
df$Fact10<-ifelse(df$FACT10==2,1,0)#FACT10 # opinion 
df$Fact11<-ifelse(df$FACT11==2,1,0)#FACT11 # opinion 
df$Fact12<-ifelse(df$FACT12==1,1,0)#FACT12 # Fact 

df$Fact.total<-
  df$Fact1+df$Fact2+df$Fact3+df$Fact4+df$Fact5+df$Fact6+
  df$Fact7+df$Fact8+df$Fact9+df$Fact10+df$Fact11+df$Fact12

df$partisan<-ifelse(df$demdum==1,"Democrat",ifelse(df$repdum==1,"Republican","Independent"
                  ))
b<-rep("lightblue",times=13)
g<-rep("lightgreen",times=11)
r<-rep("pink",times=12)
col<-c(b,g,r)
df%>%
  filter(!partisan=="NA")%>%
  ggplot(aes(x=Fact.total))+
  geom_bar(aes(y=after_stat(prop)),fill=col)+
  facet_wrap(vars(partisan))+
  theme_classic()



###############
# table 1
describe(df$Y)


###############
#table 2
df%>%
  filter(repdum==1)%>%
  dplyr::select(accurate_response,partisan_error,unbiased_error)%>%
  describe(.)


###############
#table 3

df%>%
  filter(demdum==1)%>%
  dplyr::select(accurate_response,partisan_error,unbiased_error)%>%
  describe(.)


###############
#table4
df%>%
  filter(demdum==0 & repdum==0)%>%
  dplyr::select(accurate_response,partisan_error,unbiased_error)%>%
  describe(.)


###############
#table5

df%>%
  dplyr::select(accurate_response,partisan_error,unbiased_error)%>%
  group_by(accurate_response,partisan_error,unbiased_error)%>%
  count(.)%>%
  print(n=22)



###############
#table6
df1<-df[complete.cases(df$demdum,df$repdum),]
df1%>%
  dplyr::select(age,polknow,currevent,analytical,educ,dem_bias,rep_bias)%>%
  describe(.)


###############
#table7
df%>%
  filter(repdum==1)%>%
  dplyr::select(age,polknow,currevent,analytical,educ,dem_bias,rep_bias)%>%
  describe(.)


###############
#table8
df%>%
  filter(demdum==1)%>%
  dplyr::select(age,polknow,currevent,analytical,educ,dem_bias,rep_bias)%>%
  describe(.)


###############
#table9
df%>%
  filter(demdum==0 & repdum==0)%>%
  dplyr::select(age,polknow,currevent,analytical,educ,dem_bias,rep_bias)%>%
  describe(.)


###############
# Table 10
stargazer(step1,type="text",digits=2)


###############
# Table 11
stargazer(step2,type="text",digits=2)


###############
# Table 12
L0<-multinom(Y ~ 1,Hess=TRUE, data = df)
LM<-multinom(Y ~ polknow+analytical+educ+currevent+age+repdum+demdum,Hess=TRUE, data = df)



L.full <- logLik(LM)
L.base <- logLik(L0)
n<-length(df1$caseid)
G2 <- -2 * (L.base - L.full)
r1<- as.numeric(1 - exp(-G2/n))
z1 <- summary(step1)$coefficients/summary(step1)$standard.errors
p1 <- (1 - pnorm(abs(z1), 0, 1)) * 2
coefs1<-exp(coef(summary(step1)))
stargazer(step1,type="text", 
          p=list(p1),
          coef=list(exp(coef(summary(step1)))),
          covariate.labels= c("Civics knowledge","Cognitive ability","Education","Current events knowledge","Age","Rep","Dem"),
          add.lines = list(c("n",nrow(step1$residuals),nrow(step1$residuals)),
                           c("Pseudo R^{2}",round(r1,2),round(r1,2))))



###############
# Table 13
LM1<-multinom(Y ~ dem_bias+rep_bias+polknow+analytical+educ+currevent+age+repdum+demdum, data = df,Hess = TRUE)
L.full.1 <- logLik(LM1)
G2_1 <- -2 * (L.base - L.full.1)
r<- as.numeric(1 - exp(-G2_1/n))
z <- summary(LM1)$coefficients/summary(LM1)$standard.errors
p <- (1 - pnorm(abs(z), 0, 1)) * 2
coefs<-exp(coef(summary(LM1)))
stargazer(step2,type="text",
          p=list(p),
          coef=list(exp(coef(summary(step2)))),
  covariate.labels= c("Pro-Dem affect pol.","Pro-Rep affect pol.","Civics knowledge","Cognitive ability","Education","Current events knowledge","Age","Rep","Dem"),
  add.lines = list(c("n",nrow(step2$residuals),nrow(step2$residuals)),
              c("Pseudo R^{2}",round(r,2),round(r,2))))
          
